Bistable collective behavior of polymers tethered in a nanopore 
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Polymer-coated pores play a crucial role in nucleo-cytoplasmic transport and in a number of 
biomimetic and nanoteclmological applications. Here we present Monte Carlo and Density Func- 
tional Theory approaches to identify different collective phases of end-grafted polymers in a nanopore 
and to study their relative stability as a function of intermolecular interactions. Over a range of 
system parameters that is relevant for nuclear pore complexes, we observe two distinct phases: one 
with the bulk of the polymers condensed at the wall of the pore, and the other with the polymers 
condensed along its central axis. The relative stability of these two phases depends on the interpoly- 
mer interactions. The existence the two phases suggests a mechanism in which marginal changes 
in these interactions, possibly induced by nuclear transport receptors, cause the pore to transform 
between open and closed configurations, which will influence transport through the pore. 

PACS numbers: 87.15.A-, 82.35.Gh, 87.16.Wd, 87.85.Qr 



I. INTRODUCTION 

Physical modeling is a powerful tool to interpret the 
complexity arising from multiple interacting components 
in a biological system. One such system is the nuclear 
pore complex (NPC). This structure mediates all trans- 
port between the cell cytoplasm and the nucleus, and 
its operation is thought to depend on the properties of 
natively unfolded proteins, nucleoporins, that are end- 
grafted inside a ~50 nm wide channel [IHI]. Molecules 
larger than ^6 nm can only pass through this nanopore 
if they are bound to nuclear transport receptors, which 
are known to interact with the nucleoporins. Though 
the composition, hydrodynamic size, and nanomechani- 
cal properties of single nucleoporins are known in great 
detail @H§], their collective behavior in the NPC is still 
heavily disputed. This behavior has been studied using a 
variety of models, with nucleoporins in a one-dimensional 
geometry [7 , grafted to a planar surface [5] and con- 
strained within a rectangular box j9] . Only very recently 
has the cylindrical geometry of the pore been taken into 
account [TUHT5] . 

More generally, there has been significant interest in 
polymer coatings in nanopores, since they can be used to 
tune the aperture of artificial and biomimetic nanopores 
and filters |14ffTT3] . Polymers end-grafted in a cylindri- 
cal pore have been studied by molecular dynamics |20) 
and Monte Carlo |21_ simulations. Depending on solvent 
quality, pore diameter and polymer structure and dimen- 
sions, there can be a rich and complex phase diagram 

E2HEI. 

In this paper, we investigate the effect of confinement 
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Figure 1. (Color online) (a) Illustration of the polymer model 
where bead n is constrained to lie on a circle such that the 
length of the bonds to beads n — 1 and n+ 1 is unchanged, (b) 
Snapshot of a Monte Carlo simulation with 40 non-interacting 
polymers of contour length 100 nm, bead diameter d = 1 nm 
and bond length b = 1 nm, end-grafted to the thick dots in a 
cylindrical pore of 25 nm radius. 



on the possible conformations of polymers within a cylin- 
drical channel. We have performed Monte Carlo simula- 
tions on a coarse-grained model of polymers end-grafted 
within a cylinder, and have furthermore developed an ap- 
proach using density functional theory (DFT) to study 
the relative stability of competing morphologies. The 
DFT free energy is based on a reference case of a tethered 
freely jointed chain of point-like beads, a novel choice in 
this context, and its mean field variational optimization 
has been carried out using a highly efficient numerical 
procedure. The resulting well-founded free energy esti- 
mates enable us to construct a phase diagram indicating 
the relative stability of different polymer configurations 
and to speculate about transitions that might be induced 
between them. 
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II. MODELS OF NUCLEOPORIN BEHAVIOR 

Nucleoporins in the NPC channel have been shown to 
separate into two distinct categories: those that form 
short globular conformations, and longer polymers that 
tend to be able to extend further away from their teth- 
ering point at the NPC rim [S]. As cylindrical confine- 
ment will mostly affect those nucleoporins with contour 
length much greater than the pore radius, we focus on 
this latter category: polymers with 100 nm contour 
length, end-grafted on a ring around the inner wall of 
a long cylinder with a radius R = 25 nm (Fig. [TJ. We 
model the polymers as freely jointed chains of beads of 
diameter d = 1 nm, with a segment length 6 = 1 nm. 
This implies a persistence length of 0.5 nm, in agreement 
with single-molecule pulling experiments on nucleoporin 
cNupl53 |6J. Excluded- volume effects are modeled as a 
hard-sphere repulsion between the beads. They are sup- 
plemented by longer-ranged attractive interactions be- 
tween the polymers, consistent with those that appear to 
operate in the NPC [IH3] • The combined bead-bead pair 
potential therefore takes the form 



W [-eexp[-(|r|-d)/A] \r\>d, 

where r is the vector connecting the centers of the beads, 
and A and e are range and strength parameters, respec- 
tively. A variety of interaction mechanisms might be rep- 
resented by an appropriate choice of d, A and e in Eq. (JT|) , 
including, for example, the hydrophobic interaction that 
is thought to play a significant role in these systems |25| . 

As a first approach to determine typical polymer con- 
figurations, we study the system by Monte Carlo (MC) 
simulations. We employ a straightforward Metropolis al- 
gorithm, in which single beads attempt moves on a cir- 
cular path defined by the constant distance b to their 
nearest neighbors, whilst remaining restricted to the in- 
ner volume of the cylinder as illustrated in Fig. [T] The 
first bead for each polymer, located on the cylinder inner 
surface, is fixed in position, whilst the last bead is free to 
move on the surface of the sphere of radius b centered on 
the penultimate bead. The restriction of constant seg- 
ment length simplifies the description of each MC move, 
but can slow down the exploration of configuration space. 
Nevertheless, all polymer conformations are accessible 
from one another. We perform simulations of 250000 
attempted moves per bead. Relaxation to equilibrium is 
confirmed by noting convergence of the system energy, 
such that we use the second half of each simulation to 
generate mean bead profiles. The simulations indicate 
an interesting range of behavior, as can be seen in Fig. [2] 
for 40 polymers, each of length 100 beads, tethered uni- 
formly around a ring. Different dominant configurations 
may be observed depending on the parameters chosen, 
though they can be roughly divided into two categories: 
conformations in which the density is peaked in the cen- 
ter, and those in which it is peaked closer to the wall. 



(a) 




Figure 2. (Color online) Snapshots taken from converged 
Monte Carlo simulations as in Fig. [TJ showing different re- 
sults for polymers that are subject to the same excluded vol- 
ume and attractive interactions, as defined by Eq. [I] with 
e = 0.1 ksT and A = 1.0 nm. (a): Condensation in the center, 
(b) and (c): Different numbers of clumps can be found when 
the polymers condense closer to the wall. 



Sometimes profiles from both categories emerge even for 
the same parameter choice, depending on the starting 
configuration. This suggests that there exist thermody- 
namically stable and metastable states for a particular 
parameter set. 

In order to categorize these phases more fully, and in 
particular to investigate the relative stability of wall and 
central configurations, a variational mean field density 
functional theory (DFT) of a many-polymer system has 
been developed. DFT provides a natural framework for 
evaluating the free energy of large numbers of interacting 
particles [55] , though certain modifications are necessary 
for it to be suitable as a model of a system of polymers. 
It shares many of the features of a self-consistent statis- 
tical field theory of polymers [37], and has the capacity 
to include finite-range interactions, as expressed in the 
pair-potential <p(r). The polymer entropy is estimated 
by solving an equivalent Brownian motion problem, in 
contrast to other approaches that are also based on the 
minimization of a free energy functional but estimate the 
entropy by numerically generating a large set of sample 
configurations |24| . Our model can be implemented using 
an efficient numerical algorithm that runs on a standard 
desktop PC. It can be used to explore the equilibrium 
behaviour of the system. Dynamical versions of DFT 
have been developed to include relaxational phenomena 
[281 125] , and our model has the potential to be extended 
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in this direction. 

The DFT model applied to a single chain of TV beads 
employs the Hamiltonian 



N-l 



N N 



i=0 j=0 ijtj 

where r; is the position of bead i, h is a function that 
constrains each segment of the polymer chain to take 
a fixed length b, and the potential 4> acts between all 
pairs of beads. In the variational mean field approach, 
we introduce an additional potential V such that H 
if + H 1 with 



N-l 



A'* 



H =^2h(r i+1 ,T i )+J2V(r i ) (3) 

i=0 i = l 

N N N 



H 



3=0 i^j 



i=l 



Hq describes a chain of TV point-like beads interacting 
with an external potential V, constrained by a require- 
ment of constant segment length, and with the first bead 
coupled to a static tethering point at ro. This system 
has thermodynamic properties embodied in the free en- 
ergy g7] 



F = -ln(y G(r ,r,JV;M)dr) , 



(5) 



in units of fc^T, where G(i"o, r, TV; [w]) is the propagator 
of a Brownian motion with contour length Nb from tether 
point r to point r inside the cylinder, evolving under the 
influence of a dimensionless external potential w(r) — 
V(r)/kBT acting as a sink; namely the Green's function 
solution [30] to the diffusive problem described by 

dG(ro,r,s;[ W ]) = ( b *_ w{r) \ G (g) 



d.s 



with initial condition G = S(r — Tq) at s — 0, and bound- 
ary conditions G = at the cylinder wall and zero radial 
gradient at the center. p(r) is the single bead distribu- 
tion function evaluated for the Hamiltonian Hq; it is a 
functional of w, and is given in terms of G as [5U] 



p(r) = 



ffds /dr'G(r , r, N - s; [w])G(r, r', s; [w]) 
fdr'G(r ,r',7V; [w]) 



(7) 



It is important to recognise that the mean field in DFT 
is a single-bead potential that is introduced to emulate 
the real polymer self-interactions as closely as possible. 
The reference Hamiltonian H describes the behavior of 
TV freely jointed beads in the potential V, and the ab- 
sence of additional bead-bead interactions simplifies its 
analysis considerably. The DFT-derived bead density 
profiles represent polymer configurations adopted in re- 
sponse to a mean field potential instead of the actual 



self-interactions. This can be a reasonable approxima- 
tion if the mean field is optimized, as can be recognized 
through use of the Bogoliubov inequality. Within such an 
approach, the best description of the interacting polymer 
system can be obtained by minimizing the free energy 



F ml = F - p(r)w{r)dr+F h 



p(r)p(r')u(r— r')drdr' 
(8) 

with respect to the mean field w, bearing in mind that 
p depends on w. Both F m { and it, the attractive part 
of the bead pair potential <j>, are expressed in units of 
ksT. The contribution to the free energy arising from 
the attractive term is derived on the basis of a random 
phase approximation. The contribution from the repul- 
sive interactions is modeled by a hard chain free energy 
-Fhc described in the local density approximation, and in 
units of fc R T. by [3Tj 152]: 



^hc = / P(r) 



1 - 



4?y(r) - 3?7(r) 2 
[l~^(r)] 2 
1\ / 2-r?(r) 
N) V2[l-r,(r)]3 



(9) 



dr 



where rj(r) is a bead packing fraction given by 7rp(r)<i 3 /6. 
The origin of Eq. ([8]) is described in detail in Appendix 
A. 

The model so far has been constructed for a single 
polymer, but a system of M polymers in a pore can be 
treated by multiplying F Q by M, by interpreting p in the 
remaining terms in Eq. ^ as the superposition of bead 
density profiles of the M polymers attached to their sepa- 
rate tether points, and by regarding the entire free energy 
F m { as a functional of a mean field w that we assume, 
for simplicity, to exhibit the cylindrical symmetry of the 
pore. 

We minimize the free energy functional ([8]) using a con- 
jugate gradient method (Polak-Ribiere) |27l 155] . on the 
basis of gradient information obtained by evaluating the 
functional derivative SF^f [p] /Sp. This has proved to be 
more efficient in this context than employing the gradi- 
ent information within a steepest descent method: this is 
discussed in greater detail in Appendix A. The numerical 
scheme involves making an initial guess for u>(r), evalu- 
ating the bead density using Eq. ([7]), thereby obtaining 
the mean field free energy through Eq. Q and using the 
conjugate gradient scheme to generate a new free energy, 
and hence a new mean field. The mean field and the bead 
density are defined numerically on a grid of points within 
the cylinder, fine enough such that the spacing does not 
affect the outcome of the calculations. 



III. COMPARISON BETWEEN MC AND DFT 

A key test of the DFT model, and of the numerical 
scheme, is to compare bead density profiles with those 
obtained by MC simulation. We look first at polymers in- 
teracting only through hard sphere repulsion, represented 
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Figure 3. (Color Online) Mean bead density profiles for Monte 
Carlo (MC) simulations (a) and DFT calculations (b) for 40 
polymers, composed of 100 beads each, tethered within a 
cylinder, and interacting only through hard sphere interac- 
tions. The plots are given for a range of radial (r) and axial 
(z) coordinates, (c) and (d) show more detailed comparisons 
of the profiles as a function of, respectively, the radial direc- 
tion in the plane of the tethering ring (0 = 0), and the axial 
direction at a radius r = 17.5 nm. DFT results are given 
by the solid lines and the statistical uncertainty in the MC 
results is indicated with error bars. 



in the DFT by the contribution in Eq. Q . The compari- 
son shown in Fig. [3]shows a good agreement between the 
two schemes, indicating that the DFT approach captures 
the essence of the excluded volume behavior. 

Next, we include long range interactions and make a 
similar comparison between DFT and MC results. This 
test is more challenging since there is now competition 
between polymer attraction and repulsion, giving rise to 
quite distinct configurations, as we saw earlier for the 
MC alone. Whilst the DFT can treat cases with differ- 
ent radial density profiles, it presently does not explicitly 
allow for azimuthal clumping. Fig. |4|a-b) illustrates two 
converged density profiles obtained from the DFT ap- 
proach, arising from different choices for the initial mean 
field. They are the counterparts to the wall and central 
phases observed in Monte Carlo simulations seen in Fig. 
[2] and use the same set of interaction parameters, namely 
e = O.f ksT and A = 1.0 nm. A detailed comparison of 
the centrally peaked profiles with MC results in Figs.[4jd) 
and (f) suggests that the radial and axial spread of the 
polymers determined from each approach are consistent 
with one another. For the wall phase, the DFT and MC 
profiles do differ, as illustrated radially and axially in 
Figs. |4jc) and (e). The differences might be due to the 
angular symmetry breaking, or clumping, observed in the 
Monte Carlo simulations, but this would not be expected 
to affect the qualitative conclusions about the stability or 
metastability of the two phases that we now explore. 

The great benefit of the DFT model is that it pro- 
vides thermodynamic properties of the interacting poly- 
mers within the pore, not just the structural properties 
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Figure 4. (Color online) Converged polymer configurations 
and bead density profiles calculated for an attractive pair po- 
tential of depth e = O.lfcsT and range A = 1.0 nm in a 
pore of 25 nm radius. The data in the left and right col- 
umn correspond to initial conditions with the 40 polymers 
concentrated at the wall and at the center, respectively, (a- 
b) DFT results with parameters corresponding to the Monte 
Carlo simulations, as a function of radial (r) and axial (z) 
position, (c-d) Comparison of Monte Carlo (with error bars) 
and DFT (smooth curves) radial profiles for z = 0, the plane 
of the tethering points. Axial profiles (e) for the wall phase 
at r — 22.5 nm and (f) for the central phase at r = 2.5 nm. 



that are available using MC. It provides estimates of the 
free energy, such that it is possible to determine which 
phase, central or wall, is thermodynamically stable or 
metastable under a range of interaction conditions. 

In Fig. [5] we plot the difference in free energy be- 
tween wall and central profiles AF p , per polymer and 
in units of fc^T, against interaction range A and strength 
e. The plot illustrates the free energy difference as a 
surface extending across regions where, respectively, the 
central phase ('Center stable, wall metastable') and the 
wall phase ('Center metastable') are thermodynamically 
stable and the other phase is metastable. The region 
in the foreground ('Wall only') denotes conditions where 
only the wall phase appears to exist. Similarly, there is a 
corresponding region where only the central phase exists 
towards the back of the diagram ('Center only'). There is 
a binodal line, or phase boundary, where wall and central 
phase coexist, together with spinodals denoting the ex- 
tremes of metastability of one of the phases with respect 
to the other. When we extend the range of e to larger 
values, i.e. stronger attractions than those shown, we 
find that the phase boundary continues in such a manner 
such that a central phase is increasingly favored. 
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Figure 5. (Color online) Metastability diagram of the polymer 
phases as a function of strength e and range A of the attrac- 
tive interactions. AF P indicates the free energy difference per 
polymer between the wall and central phases, for parameters 
where they can both exist as a stable and metastable state. 
The coexistence conditions lie along the boundary between 
blue and red regions. In the green region, only one phase is 
found to be possible: the metastability of the other has been 
lost. 



The metastability diagram demonstrates that central 
polymer condensation is a natural result of attractive in- 
teractions, provided that their range and strength are 
sufficiently large. This is in accord with intuition, which 
suggests that a central phase can be stabilized by a reduc- 
tion in energy to balance the cost in entropy of extending 
the polymers away from the wall. A long range attractive 
potential will favor this by allowing polymers to interact 
with more of their counterparts across the other side of 
the pore. If the range were reduced, then the required 
strength of the attraction would have to be greater to 
produce the same effect. Repulsive interactions do not 
drive bistable phase behavior. Instead the polymer con- 
densate becomes more homogeneous as the strength of 
the repulsion is increased: such behavior is equivalent to 
that of a good solvent. 

In addition to changing the strength and range of the 
polymer interactions, we can also explore the effect of a 
change in the number of beads N in each polymer. For 
constant segment length b, this is equivalent to altering 
the length of the polymer. Reducing TV has the dual 
effect of increasing the entropic cost of stretching a poly- 
mer towards the center of the pore, and of decreasing the 
binding energy that can be experienced by each poly- 
mer. Both effects hinder central phase formation, and 
this is borne out in calculations of free energies for the 
two phases as a function of interaction strength e for a 
range of N, as shown in Fig.[6]for A = 0.5 nm. Shortening 
the length of the polymers requires a stronger attractive 
potential for the polymers to condense towards the pore 
center, as is to be expected. 
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Figure 6. Free energies of the central phases (solid) and wall 
phases (dashed) for systems with a range of bead number 
N, per polymer, for interaction range A = 0.5 nm, plotted 
as a function of interaction strength e. The phase with the 
lower free energy is thermodynamically stable. The limited 
extent of the curve of higher free energy, in some cases, is an 
illustration of the limits of the metastability of that phase, 
equivalent to spinodal behavior. As N increases, the central 
phase becomes stable over a larger range of e. 



IV. DISCUSSION AND CONCLUSIONS 

On the basis of Monte Carlo and density functional 
theory calculations, we find that polymers tethered 
around a ring inside a cylindrical geometry can exhibit 
bistable behavior, switching between a wall- and a cen- 
trally condensed phase depending on the interaction pa- 
rameters. Interestingly, this behavior is observed for a 
geometry and polymer structural properties that closely 
resemble the NPC [lH3|, for entirely realistic ranges 
(A < 1 nm) and strengths (e < 1 fc^T) of intermolec- 
ular interactions [25] within the NPC. The existence of a 
central polymer condensate is reminiscent of the central 
'plug' or 'transporter' structure of nucleoporins that has 
been observed by cryo-electron microscopy of the nuclear 
pore complex (NPC) [H[S], though our model would need 
to be developed to treat a more realistic geometry if it 
were to be considered a proper representation. Never- 
theless, our thermodynamic model suggests that a phase 
transition from the central to the wall phase occurs when 
the effective strength or range of the attraction between 
polymer segments is decreased. Equivalently, the transi- 
tion might be induced by increasing the prevailing tem- 
perature. If we were to apply our model to the NPC, 
this feature would be in agreement with the experimen- 
tally observed dissolution of the central plug at higher 
temperature [33]. A similar effect has been seen upon 
incubation of the NPC with nuclear transport receptors 
|35j . A possible interpretation of this effect is that the 
nuclear transport receptors are responsible for a weaken- 
ing of the nucleoporin self-attraction and that their infiu- 



6 



ence on the polymer plug provides a mechanism for the 
differential permeability of the NPC. In this scenario, a 
sufficiently high concentration of nuclear transport recep- 
tors disturbs the thermodynamic stability of the central 
structure to the extent that the polymers withdraw to- 
wards the wall, leaving a free passage to be occupied and 
traversed by receptors and receptor-bound cargos. Such 
a mechanism, if it operates, could be exploited in the 
design of artificial nanopores that might perform similar 
differential transport of a variety of molecular species. 
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Appendix A: Details of the free energy functional 

We provide here some theoretical and numerical back- 
ground, including a detailed specification of the bead 
density profile, the rationale for the variational prin- 
ciple employed, and the condition for minimizing the 
free energy using functional derivatives. We focus our 
interest on the statistical properties of a Hamiltonian 
Hq for a freely jointed polymer of N point-like beads 
at positions {r^} in a dimensionless external potential 
w, given by Hq = U ({rj) + k B TY^=i w{*i), where 
U {{^i}) = X^z M r i+i) r i) i s a function that imposes the 
constraints on bond length. The bead at ri is constrained 
with respect to the tether point at rQ. The partition func- 
tion of the system is 



N 

Zq = / \\ dr,j exp 



N 



-U/k B T -J2 w ( r i 



i=l 



(Al) 



The configuration-dependent bead density over con- 
tinuous spatial position r is defined as p(r, {r^}) = 
J2iLi — r i)' which allows us to write Z as 

Z [w] = ^[dr,] exp (-U/k B T - f p( T )w(r)dij , (A2) 

where [dr.;] represents the integration over bead positions 
and we compress the notation of p for clarity. 

The partition function is clearly a functional of the 



external potential w and its functional derivative is: 
SZ [w] _ 
Sw(y) 
1 



lim 



dr,-] exp 



U 



p(r)(w(r)+eS(r — y))dr 



[dr, ] exp 



U 



p(t)w(t) dr 



= - j [dr i ]p(y)exp (-U/k B T- j p(v)w(r) dr^j , 

(A3) 

such that we can define p(y) = (p(y)) to be the mean 
bead density for the system: 

P(y) = yJ [dr i ]/5(y)exp (-U/k B T- j p(r)w(r)dr 



1 6Zq [w] 
Z [w) 6w(y) 



S In Zq [w] 
Sw(y) 



(A4) 

which demonstrates that p is a functional of w. The 
explicit functional dependence is given in Eq. Q. 

Now we discuss a self-interacting polymer described by 
the Hamiltonian 

1 N N 

incorporating a potential <j) acting between all bead pairs. 
We write H = H Q + H x with 

N N N 
j = l i^tj i=l 

and employ the Bogoliubov inequality F < Fq + (Hi), 
where F is the free energy of the system described by 
Hamiltonian H . and Fq = — In Zq is the free energy of 
the reference system described by Hq, in units of k B T, 
and is given by Eq. |5| [27]. As before, the brackets 
denote an average over the ensemble associated with Hq, 
and we write 



z\> j ^ / ^ y ) w ( y ) d y ex p ("fcgT 



p(yMy)dy. 



p(y)w(y) dr 



(A7) 

In a similar fashion the mean of the pairwise terms is 

INN \ 
\i=l i#j / 

Y J [dr»] j P2(x, y, {r fc })0(x-y)dxdy 

-Jp(r)w(r)drj =Jp 2 (x,y)<j)(x - y)dxdy, 

(A8) 



x exp 



U 
kuf 
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where p 2 (x, y, {r fc }) = J^ytffc ~ r i) S (y ~ r j) is thc 
two-point configuration-dependent bead distribution and 
P2(x,y) = (/5 2 (x,y, {r fc })) is its mean in the H ensem- 
ble. For simplicity, we take a random phase approxima- 
tion and represent pa(x,y) by p(x)p(y). Thus the free 
energy of the self-interacting polymer is bounded by the 
inequality 



F< - In Z [w]-j p(y)w(y)dy+- / p(x)p(y)0(x-y)dxdy, 

where <f> is the pair potential divided by fc^T, which can 
be separated into an attractive part u(x — y) and a repul- 
sive part. The latter's contribution to the right hand side 
may be represented by a functional -Fhc[p] (Eq. ([9|) that 
has been found to capture the thermodynamic properties 
of hard chains: freely jointed polymers of finite size hard 
spheres [311 I3"2] . The free energy then satisfies 

F < F mf [w] = - In Z [w] + F hc [p] - /p(y; [w])w(y) dy 



1 



p(x; H)p(y; [w])u(x - y) dxdy, 



(A10) 

where the functional dependence of p on the mean field 
w is explicitly noted. 

The best estimate of F is identified by functional 
minimization of the mean field formulation -Finf[w] over 
all possible w, to be achieved by setting the functional 
derivative SF m f/Sw(r) to zero. Several contributions to 
the derivative arise. We already have S\n Z /Sw(r) — 
— p(r; [w]) from Eq. (A4|, and furthermore, regarding 
-Fhc as a functional of either w or p, 

^hcM = f SF hc [p]dp(y; [w]) = t ' 6 'pjy; [w]) 
8w(r) J 5p(y) Sw(r) J c Sw(r) 

(All) 

where /ih c represents the functional derivative of F^ c with 
respect to p, together with 

^4)/p(y; M)^(y)dy= P (r-, N)+/ <5/ ^r-(y) d y 

(A12) 

and 



5w(r) 



p(x; [w])p(y; [w])u(x - y)dxdy 



= 2/M^Mp(x;MHx-y)dxd y) 

giving the minimization condition as: 
SF mt [w] [8p(y;[w 



(A13) 



Sw(r) J Sw(r) 
5p(y; [w]) 



w (y ) dy + / p hc (y ) s _ f _^ dy 



Sw(r) 



+ 



Since 



6w(r) 



SF mf [w] 
5w(r) 



p(x; [w])u(x. - y) dxdy = 0. 



SF m{ [p}Sp(y; [w]) 
<5p(y) Sw(r) 



dy, 



(A14) 
(A15) 



this is equivalent to the condition 



w{y) + Phc(y) + / P(x)u(x - y)dx = 0, 



SF ml [p] 

S P (y) 

(A16) 

which can be regarded as a requirement that the optimal 
mean field acting on each bead is a suitable embodiment 
of the pair wise interactions. 

We now discuss algorithms to determine thc optimal 
mean field w and associated density p that satisfy this 
condition. A steepest descent method could be employed 
such that p is updated incrementally and repeatedly in a 
direction down the local slope of the F m [[p] surface. This 
can be regarded as an evolution of p in a fictitious time 
t according to 



dp(y,t) SF ml [p] 



Of 



s P (y,ty 



(A17) 



until convergence at a time-independent density profile 
where the left hand side vanishes. But the right hand 
side of this equation requires w as a functional of p. which 
is not readily available. It is easier to determine p for a 
given w, through Eq. ([7]), and we might therefore con- 
sider a scheme 



dw(r,t) _ SF mi [w] 
dt 6w(r, t) ' 



(A18) 



but the problem here is that the right hand side, given 
by Eq. (Al5l, requires a specification of the functional 



derivative <5p(y; [w])/5w(r). Instead, we employ the fol- 
lowing arguments to formulate a third scheme. From thc 
definition of p, we can write 



Spjy, [w]) 
Sw(r) 

1 



Sw(r) 



1 SZ [w] 
Zq[w] 6w(y) 



5Z [w] SZ [w] 



1 



6 2 Z {i 



{Z [w}) 2 <M r ) Sw(y) Z [w] 5w(r)5w(y) 



, (A19) 



and in view of Eqs. ( |A3| and (A4) this becomes 
Using the definitions of p and p 2 this gives 



(A20) 



P( r )/>(y)-P2(r,y)-(^£(r - r 4 )<5(y - r,-) 



$p(y; M) 

5w(r) 

(A21) 

The random phase approximation that we have employed 
asserts that pa(r, y) = p(r)p(y), so we have 



My; M) 

5w(r) 



^2S(r-r l )S(y-r l 



(A22) 



and we reach the important conclusion that, if viewed as 
a matrix, to this level of approximation the off-diagonal 



8 



elements of Sp(y; [w])/6w(r) are zero, and the diagonal 
elements are never positive. 

Now we define a functional F[w] that satisfies 



SF[w] 
Sw(r) 



SF m i[p\ 
' Sp(r) 



(A23) 



and compare the variational properties of F with those of 



F m f . In view of Eq. ( A15 ), the mean field at a stationary 
point of F, where 5F[w\/8w(r) = 0, will minimize F m {. 
Next we establish that an increment in w in the direc- 
tion of steepest descent of the functional F gives rise to 
a decrease in the value of the functional F-n& . We show 
this by evaluating the analog in function space of the 
dot product betw een t wo gra dient vectors. We calculate, 
using Eqs. (A22l and (A23l, the projection of the func- 
tional derivative of -FmflwJ along the direction in w space 
of the functional derivative of J^w], namely the integral: 



/ 



SF[w] SF ml [w] 
5w(r) Sw(r) 



dr = — 



*SF mf [p] 5F mf [p] 6p(y) 
M r ) Sp(y) 8w(r) 



I N 

E 



5F mf [p] 5F mi [p} 
5p{vi) 5p{vi) 



dydr 



(A24) 



so that J[SF[w]/6w{r)][6F m {[w]/6w(r)]dr > to this 
level of approximation. The functional F has a mini- 
mum at the mean field w that minimizes F^, and an 
infinitesimal move along a path of steepest descent of the 
F surface also takes us downhill on the F m f surface. 

All this provides support for an algorithm for identi- 
fying the optimal mean field and free energy based on 
solving the equation 



dw(r, t) 



SF[w] _ SF mi [p] 



(A25) 



dt 8w(r, t) Sp(r, t) ' 

which is free of the above-mentioned problems of schemes 



(A17I and (A18) 



Putting this scheme into practice, and using Eq. 
( A16 1, we discretize the spatial variable and the time 



and employ a forward Euler numerical scheme with up- 
date rule 



w n+i (y k ) = w n (y k ) + At -w n (y k ) + M £ c (y fc ) 



53p n (x,)«(x i -y jfc )Ax 



(A26) 



where n labels discrete time and subscripts j and k label 
discrete spatial points, and Ax is the volume element. 
At each iteration starting with a mean field w n (yk), we 
use Eqs. (6) and (7) to generate the reference system 
bead density profile p n (y ) that is associated with this 



choice. Through Eq. ( A26 1 with a given timestep At this 



gives us a new mean field w n+1 (yk) and the process is 
repeated until the change in mean field falls below a cho- 
sen tolerance. The converged field gives a minimized free 
energy F m { which provides an upper limit to the actual 
free energy F of the self-interacting polymer system. 

In fact, it proves to be more efficient to conduct the 
minimization of -F m f[w] using a modified conjugate gra- 
dient scheme. It is well known that such a scheme nor- 
mally takes the form of repeated line minimization of the 
functional along directions chosen in w space that are 
selected according to the local gradient 8F m f/8w(r) and 
the direction of the preceding line search. We employ the 
Polak-Ribiere version of this scheme but our modification 
is to select directions based on the functional derivative 
SF/Sw(r) instead. This is in the same spirit as the use 
of Eq. (A25 1 instead of Eq. (A18). The minimization of 



F m { is performed numerically by stepping along the cho- 
sen direction in a space spanned by the discrete w n (y k ) 
until we encounter a change in sign of the difference AF m f 
with respect to the previous value. A new search direc- 
tion is then established and the process repeated. The 
scheme appears to be numerically robust in practice, an 
indication that our consideration of the properties of the 
functionals _Fi n f and F is sound. 
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